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1. CURRENT TECHNICAL OBJECTIVES 




1. Optimal Utilization of Laser and VLBI Observations for Reference Frames 
for Geodynamics 

2. Utilization of Range Difference Observations in Geodynamics 

3. Estimation Techniques in Crustal Deformation Analysis 
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2 . 


ACTIVITIES 


2.1 Earth Rotation Parameter Determination from Different 
Space Geodetic Systems 


Introduction 


Since the last report on this study, much work nas been done on 
developin'^ and/or -testing; software to allow simulation and adjustment of Lunar 
Laser Ranging(LLR), Satellite Laser Ranging (SLR), and Very Long Baseline 
Interferometry (VLBI) data. This progress on readying the software for each 
jystem is described separately below, adong with some comments on the 
adjustment program to be used and other work to be done. 


SLR 


The NASA GEODYN program (version 8210.1) has successfully been used to 
simulate SLR data. Given approximate station coordinates, observational 
accuracies, observation cutoff elevation angle, time period, a satellite state 
vector, and earth rotation information, there is no problem in simulating any 
amount of data. Although the effects due to other systematic error sources 
existing in SLR could be easily considered, at least initially a simple point-mass 
earth and satellite model will be used. The operation of this specific option 
has also been verified. 

The problem noted in the Iasi report that this program was failing due to 
0C4 errors has been avoided by not printing the solar flux tables (with a 

"FLUX 2" card). This temporary fix to the problem •was provided by Dr. Pet» r 

Dunn of EG&G and was in turn reported to Ms. Barbara Putney of NASA/GSFC. 


LLl 


To simulate LLR data, the same basic assumption mentioned above will be 
made, namely, that the satellite (in this case the moon) will be considered as a 
simple point mass body revolving around a point mass earth. The situation is 
then identical to that of SLR, except that an approximation of the moon's orbit 
will be used rather than that of Lageos. As before, the main reason this 
simplifying assumption is made is to allow the simulation to be as 
straightforward as possible, by ignoring the systematic effects which are not 
considered likely to propagate strongly into the earth rotation parameter (ERP) 
estimates. In addition, there is no software currently available here to 
consider some of the systematic effects occurring in LLR, such as the lunar 
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physical librations. (No further indication has yet been received as to when 
GEODYN will be changed to allow real LLR observations; however, as mentioned 
in the last report, this is apparently still a lon<<-term ^oal of GSFC.) 


VLBI 


The simulation of VLBI data is being pursued along the following lines: 

(1) IRIS (the International Radio Interferometric Surveying Network) VLBI 
schedules have been obtained from thn National Geodetic Survey (NGS) (via 
Ross McKay on June 18, 1984). These schedules are in the VLBI Mark III 
System SKED program format as described in [Vandenburg and Schaffer, 1983] 
(also received from Mr. McKay on August 22, 1984). These schedules are, the 
actual schedules now in use by the five IRIS observatories. Since in the next 
few years substantial changes in the number of such o': -ervatories 
(permanently dedicated to ERP observations) and in the observing schedules 
seem 'unlikely, these schedules will be used to define all the VLBI observations 
to be simulated. The only change will be to consider various observation 
session intervals, rather than only every five days. And if substantial 
changes should occur in the schedule or number of observatories, the new 
schedules could be easily obtained from NGS and used instead. 

(2) Since GEODYN has never been properly set up to simulate VLBI 
observations (as it has for SLR), an intermediate program is being written 
("SKEDVIP") to read the schedules and provide information to another program 
for data simulation. A modified versioi. (now also being made) of Yehuda 
Bock’s "VIP" program (Bock, 1980] will be used to generate? these simulated 
observations in a GEODYN-compatible format. 

(3) As reported in the last semiannual report, the current version of 
GEODYN (8210.1) will not process VLBI observations since the changed source 
regarding VLBI which was to have been installed in that version was left out. 
Discussions with Dr. Dunn and Ms. Putney in August have indicated that these 
changes will also not now be added in the future but may be incorporated 
eventually in the new GEODYN II program. (GEODYN itself has been 
permanently frozen at the 8210.1 level.) 

Rather than wait an extremely long time for GEODYN II to become available 
with those changes, the changed source itself was requested from their author, 
N. Zelensky of EG&G so they could be incorporated here into a special version 
of GEODYN. One attempt at sending these changes failed in late September due 
to a tape mixup at GSFC, and another tape is awaited. When it arrives, the 
special version of GEODYN will be easily created and tested with a VLBI test 
run previously sent by Mr. Zelensky. 

(4) In any case, once the program is operating here with VLBI data, as 
with SLR and LLR it can provide individual solutions directly or pass the 
normal equations to the SOLVE program for individual or combination solutions. 
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The SOLVE Program 


The SOLVE program (version 8212.0), which will be used at least to 
provide the combination SLR-LLR-VLBI solutions, has now been tested with 
simulated SLR data. Although apparently working, there are still some 
problems occurring, wiich are still being checked to verify the correct 
operation of the program. 


Summary and Future Work 

The readying and testing of the software to simulate and adjust the three 
d'lta types is now well underway and should be completed by the time of the 
next report. 

• 

Final decisions on the stations participating, observational accuracies, and 
reference frame differences to be used should also be made soon. The main 
problem to be considered next is the simulation of the ERP information itself, 
including over what time period they are to be estimated. Software to combine 
the estimates of the individual systems using a weighted averaging method (a 
la BIH) also will be obtained or developed. 


References 


Bock, Yehuda (1980) "A VLBI Variance-Covariance Analysis Interactive Computer 
Program,” Ohio State Univ. Dept, of Geodetic Science and Surveying 
Rep. 298, Columbus. 

Vandenburg, N.R. and D.B. Schedfer (1983) "Standard Schedule File Format - 
Mark III Software Documentation,' NAJ'^ Goddard Space Flight Center, 
Greenb-jlt, Md. 



2.2 


UtiUzation of Range-Difference Observations 
in Geodynamics 


Introduction 

The research reported in the last seminannual report in reference to the 
utilization of simultaneous laser range-differences for the determination of 
baseline variations evidently demonstrated that a 30 m disagreement between 
the GE0SPP81 and GEODYN baseline estimate came from the presence of 
blunders in the original laser ranging data set. Therefore, at that point no 
conlcusions could be drawn regarding the 'performance of the simultaneous 
laser range-difference method as applied in the determination of baseline 
variations. Naturally, rejection of all the observations containing blunders 
should follow in this investigation. 


uata Processing 

The rejection of the observations containing blunders was accomplished 
through the so-called "data-snooping" procedure. The performance of this 
procedure for editing laser range observations is well established (see 12th 
Seminannual Status Report) because rejection of less than 10% of the 
observations out of the original data set first randomizes the residuals and, 
second, reduces their magnitude between -90 and 90 cm. Before applying this 
procedure the residuals ranged between -300 and +300 m. Based on the fact 
that 1% of the residuals had magnitude greater than 75 cm and the remaining 
ones ranged between -40 and 40 cm (see Fig. 1), it was decided to edit out the 
observations with residuals greater than 75 cm. A sample plot of the resulting 
residuals is shown in Fig. 2. This plot confirms that the observations which 
survived in the above procedure form a "clean" (i.e., no presence of blunders) 
laser-range data set. 

Based on this clean data set, the next step is to create the 
quasi-simultaneous ranges by first determining which of the two stations — 7114 
or 7115 — has the denser data distribution for each pass. Once this is 
determined, the data of the station with the most observations are fed into an 
interpolator, and ranges are obtained for each of the data points in the 
alternate station’s record. These ranges, together with the corresponding ones 
of the alternate station, are subtracted to produce simultaneous 
range-difference observables which are part of the input for the final 
adjustment program (i.e., GE0SPP81). It is clear, therefore, that the choice of 
the interpolation method is very critical, because choice of the wrong 
interpolator will annihilate the precision of the third-generation lasers. Among 
the interpolation methods, for reasons well explained in [Pavlis, 1982], the cubic 
spline and Chebychev interpolators were chosen. Through these two 
Interpolators, the quasi-observables (simultaneous range-differences) were 
generated for two sample passes (i.e., pass *11 and pass *12) and together 
with 


— Lageos initial state vectors for the starting epoch at each pass 
— short planetary ephemeris file for the period span of the observations 
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Fig. 1 Residuals of Chebychev interpolation (after data snooping) 
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— coordinates of the pv.ie and variations in UTl at the epochs of the 
observations 

were fed into the GEOSPP81 adjustment program. The coordinate system in 
these solutions was realized through the following standard deviation valuee. 

— station 7114, fixed (i.e., <rx = vy = <rz = 0.0001 m) 

— station 7115, ax = vy = ffz = 5.0 m 

— initial state vectors: position, <rx ~ vy z az ~ 10.0 m 

velocity, vx : <ry zoz - 10 cm/s 

The orbit residuals as computed using the quasi-observables (i.e., 
simultaneous range-differences) which were generated through cubic spline and 
Chebychev interpolators have an rms of 22 cm and 11 cm respectively, and the 
plots of these residuals are depicted in Figs.- 3 and 4. Inspecting Pig. 2 the 
existence of a high correlation between orbit residuals (GEOSPP81/CUBIC 
SPLINE) and gaps in successive data points can be seen. This, of course, is 
very likely to happen because the cubic spline interpolator is a piece-wise 
interpolator and therefore very sensitive to the successive data gaps. 

To confirm this sensitivity, the differences between cubic spline 
quasi-observables (simultaneous range differences) and Chebychev 

quasi-observables were computed and plotted (see Pig. 5). These differences 
have exactly the same behavior as the orbit residuals (GEOSPP81/CUBIC 
SPLINE) (compare Figs. 3 and 5), and this assures us that the high correlation 
of the orbit residuals (GEOSPP8 1/CUBIC SPLINE) and the successive data gaps 
comes merely from the cubic spline interpolator in view of the fact that the 
Chebychev interpolator is a global, gap insensitive interp>olator. Hence, cubic 
spline interpolation should not be used when successive data gaps are greater 
than 10 s. This is exactly the situation with the data set under question. 
Therefore, the Chebychev interpolator was deployed to generate the 
quasi-observables (simultaneous-range differences). 

Next, using these observables and in response to D. Christodoulidis’ 
request in July, 1984 (private communication), it was decided to compute 
through the GEOSPP81 program the 7114-7115 baseline for the periods shown in 
Table 1. Th-.je baseline estimates can then be compared with the 
corresponding ones as determined by the GEODYN program at GSFC 
Geodynamics Branch and which were k'ndly supplied to us by D. 
Christodoulidis (private communication, August, 1984). 

Before estimating the 7114-7115 baseline, its relative orientation with the 
satellite passes for the periods shown in Table 1 should have been examined 
(see Pig. 6). Fortunately enough, 95% of the passes for each of these periods 
are parallel to the estimated baseline, assuring, therefore, its reliable 
determination. Inspecting Table 1 it is clear that both methods (i.if., 

simultaneous range-difference mode and range mode) perform about the samc^ 
with a bias of about 17 cm. It is very interesting to see that the rate of 
change of the baseline in successive periods is about the same in both 
methods. But does it really mean anything? 
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Fig. 3 Orbit residuals after data snooping and editing (750 cm residuals) 
(GEOSPP81/CUBIC SPLINE) 
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Concluaiona 


It ia evident from the above diacuaaion that the GEOSPP81 baaeline 
eatimate ia compatible with the GEODYN one, but the question of the bias 
exiatence ahould be anawered aa well as the performance of the simultaneous 
range difference method in determining real (due to crustal motions) baseline 
change, over extended periods of time, for instance one or two* years. This is 
exactly the direction which this investigation will follow. 


References 


Pavlis, E.C^ (1982) "On the Geodetic Applications of Simultaneous 
Range-Differencing to Lageos," Dept, of Geodetic Science and 
Surveying Rep. 338, Ohio State Univ., Columbus. 


2.3 An Algorithm for Crustal Deformation Analysis 


Introduction 


The complexity and unknown nature of the laws of crustal deformations 
have forced geodesists to approach geotectonics problems only by providing 
the necessary measurements to related disciplines at the beginning. Although 
the physical interpretation of the underlying mechanism of deformation still 
comes from other disciplines, a more detailed analysis of deformations becomes 
a necessity. This is caused mainly by the departures of observed quantities 
from the expected deformations along plate boundaries predicted by the 
long-term averages of crustal deformations. This led geodesists to modify their 
simplistic approach of providing only the changes in repeatedly observed 
geodetic entities to more complicated ones. 

Currently the problem is to find models of the deformation process which 
fit well to the measurements and explain the discrepancies between small scale 
and global components of driving mechanisms. Some of the difficulties in the 
analysis in both scales may arise from single modeling. Therefore, elucidation 
of deformation mechanisms is, in nature, inescapably iterative on different 
models. 

Today there exists an inherent propensity among geodesists to employ the 
sample oriented analysis of deformation measurements. However, increasing 
measurements indicate a better agreement with theoretically predicted values. 
It is, therefore, another kind of inference, namely, improving and testing global 
solutions through sample information rather than inferring them from local or 
regional studies, should also be kept in mind. 

Fig. 1 depicts an algorithm that can be used for the analysis of regional or 
local deformsition measurements. Some new techniques and testing procedures 
are proposed. 


Model Identification and Improved Estimation 

Geodetic observations are subject to error. In the first step they must be 
tested for blunders and outliers. Diagnostic checking is performed on each 
epoch of observations. Since observations are invariant quantities with 
respect to datum definitions, a proper minimum constraint solution is possible. 
Prior information on parameters or on observations are most likely to increase 
residuals and effect sample statistics. Therefore adjustments at this stage 
should always be based only on sample information. 

The second step in the analysis is the identification of the proper 
descriptive model of deformation process. The qualitative information provided 
by other disciplines and previous in situ measurement analysis play an 
important role at this stage. Nevertheless, postulating a general class of 
descriptive deformation models is mostly based on the experience and intuition 
of the analyst. Recently, different global tests have been developed for 
discriminating different models [Chrzanowsky, 1981]. In this study, a new test 
based on the measure of entropy is being investigated for discriminating 


models from the observations performed at different epochs. Several models 
are tested until an agreement is reached between model and data. 

In the next step quantitative prior information is identified and tested 
versus the sample information (Appendix B). If prior and sample information is 
found to be compatible, an improved estimation of parameters is performed. 
Different estimation techniques are being tested for improved estimation. An 
efficiency function is set up to measure the gain for different estimators 
aga..ist weighted least squares estimation (Appendix A). 


Referenc e 

Chrzanowsky, A. (1981) "A Comparison of Different Approaches into the Analysis 
of Def -mation Measurements," Proc. of XVI FIG International Congress, 
Switzerland. 
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APPENDIX A: CRITERION FOR COMPARINH ESTIMATORS AND PREDICTORS 


Comparini; two estimators (particularly unbiased versus biased), we have to 
weigh the advantages against the disadvantages of both estimators. Obviously, 
there exist so many ways to express such preferences. Each can be best for 
different purposes. Unfortunately, a universal comparison criterion does not 
exist for such purooses. As a criterion of "betterness" the scalar mean square 
error criterion is considered. The reasoi for such a choice is discussed 
through the definitions. An efficiency function based on the sbalar mean 
square error is set up. 

Definition: Mean Square Criterion (MSE). An estimator of ( is said to be 

strongly better than i^ the difference 

A(?,. ?a) = MSE(?a) - MSEC?,) 

is non-negative definite. 


Definition: Scalar^ Mean Square Error Criterion (SMSE). ( , is said to be 

SMS E-better than ?a if it has smaller scalar MSE, i.e., if 

E(?a-0T(?a-0 " E(f,-OT(?,-n 
= tr{MSE(?a) - MSEC?.)} ^ 0 

Clearly, if ?, is MSE better than ?a, then this holds for the SMSE criterion 
too. 


The problem of optimal estimates can only be attacked after the 
establishment of a generally acceptable criterion of optimality. Specific 
practical problems are sometimes o«st treated> with different criteria of 
optimality. 

A 

Let L((, O be a function which, for each true value of perrameter ( and its 
estimated value (, assigns a "loss" or "cost." The loss function cannot always 
be chosen as one which most naturally fits the physical problem. Compromises 
in the nature of the function are usually necessary, and the selection of a loss 
function in a specific case is often governed more by manipulative convenience 
than by physical eu'guments. 


Definition: A loss function < a member of class L, I ^ L if it has the following 

properties 


•L: l(q) = 0, 

<(qa) ^ <(qi) 
<(q) = <(q) 


if q = 0 

if qa * q, » 0; monotonic 
symmetric 


where q = (-(. It is generally convenient to choose the loss function as 
convex function (for minitnization purposes). A natural selection of a loss 
function from the mathematical and physical point of view is 

L(?, A) = A(?-() 

where A is a positive definite symmetric matrix of corresponding order. 


Definition: Risk Function. Given the loss L(?, A), the risk R is defined as the 

expected value of the loss function, 

R(?, A) = E(L(?, A)) = E{(?-()T A(?-n> 

Since a priori restrictions that reduce the dimensionality of the parameter 
space results in an increase in estimation efficiency, then the question is how 
important is this increase in practice. That is, we need to investigate the 
"magnitude" of the efficiency increase that results from a priori restrictions. 
In order to speak of the amount by which efficiency is increased, we must 
define some cardinal measure of efficiency. 


Definition: A measure of efficiency of an estimator is defined by the following 

scalar 

6 = (R(?i) - R(?a)) / R(?i) 

which gives the decrease of risk relative to the sample estimation. The nearer 
6 is to one, the higher is the influence of the prior information on the 
efficiency of the estimation. If 6 tends to zero, this indicates that the sample 
Information is the dominating one. If the risk is chosen as the SME, the 
efficiency function is by the following expression 

6 = {tr MSE(?i) - tr MSE(?a)} / tr MSE(?,) 
or 

6 = {SMSE(?,) - SMSE(?a)} / SMSE(?,) 


APPENDIX B: TEST FOR COMPATIBILITY OF PRIOR AND SAMPLE INFORMATION 


Once the model of deformation is identified (discriminatory analysis), the 
next step is to check the possibility that prior and sample information are in 
agreement with each other. This step is one of the deciding factors on 
whether prior information should be included in the estimation procedure. 

Let X denote the prior information on deformation parameters J-H+e where "S 
is the error vector. Then the following null hypothesis is tested: 

H : Prior and sample information are in agreement. 

Under this hypothesis if prior information and an estimate of parameters of 
WLS type are in agreement, then their difference is expected to be close to 
zero, 

6 = X-? =e-N“‘ATw-‘H 

where N = (a’^W”‘A), W the weight matrix, and u is the error vector of 
observations. The matrix of second moments of this difference is 

E(<J<5’) = Eee - 

if 

E(oi ’ ) = 

where z = ^ee~ ^ee the covariance matrix of a priori information, 

is the a priori variance of unit weight. Then, 

«Tz-% - (0, I) 

If the disturbance terms in sample observations and prior information are 
assumed to be distributed normally, i.e., 

n ' 11(0, <r» Euu), Z ~ N(0, Eee) 

then 

' N(0, I) 

It can be shown that the scalar T = can be used as a lest statistic 

since it follows a central x* distribution. T is called "compatibility statistic" 
[Theil, 1963], 

T = «Tz-% s 

= (X-T)T (Eee + s* N-‘)07-T) 

can be used which has asymptotically the same distribution. Here, s is an 
estimate of a priori variance of unit weight. 
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Theilt H. (1963) "On the Use of Incorapolete Prior Information in Regression 
Analysis," American Statistical Assoc. J. . 58, 401-414. 


2.4 Preliminary Results of the Reference Frame Comparisons in Torms 
of the Pole Coordinates from Different Techniques 


The data sets used in the present comparisons are from Project MERIT 
Circulars issued by the Joint Working Group on the Rotation of the Earth and 
from the Earth Orientation Bulletins issued by IRIS. The data involves 
astrometry (BIH, IPMS), Doppler (coded as Dop92( 1967-921), Dop67( 1970-671) and 
Dop44( 1981-441)), SLR (Lageos, coded as CSR39) and VLBI. Starlette is 
excluded from comparison because of its sparse amount of data. All the data 
sets, except Dop92 and VLBI, span the whole comparison period of September, 
}983, through June, 1984. The Dop92 data for September, 1983, February, May 
and June, 1984, are not available. The VLBI data for September, 1983, are also 
not available. 

The data are smoothed. The smoothing models of the coordinates of the 
pole are as follows: 

Model 1: 
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where 

h = 2n (MJD - 45578) / 365.25 
C = 2« (MJD - 45578) / 435 

in which MJC 45578 corresponds to September 1, 1983. 

Table 1 lists these coefficients with their standard deviations and post-fit 
standard deviations of "observations." The values of coefficients k, and k« 
representing relative polhode centers are drawn in Fig. 1. It should be noted 
that the time intervals of tabulated values vary among different techniques. 
Here the five-day values for BIH and IPMS data, the four-day values for 
Doppler data, and the roughly five-day values for SLR and VLBI data are 
used. 

Table 1 and Fig. 1 show that no large systematic difference exists in the 
pole origin. Note that the BIH coordinates and the coefficients derived from 
them have larger standard deviations. As mentioned earlier, Doppler 92 lac ^ 
four-month data, so it has a relatively large standard deviation. 

For comparison purposes, the means of pole coordinate differences between 
various techniques are computed, using the simple formulae: 


9 


n 

1 

n 

‘Via = hI (y^ - y^^) 


Table 2 includes these means of differences along with the differences 
between k|’s and between k«’s as listed in Table 1. It should be noted that in 
order to obtain simultaneity among all types of "observations" from different 
techniques the five-day pole coordinates for CSR39, Doppler 44 and VLBI are 
interpolated at the time of BIH coordinates, using a fourth-order Lagrangian 
interpolation method. It is found that Axm and Ak,, and Aym and Ak« tend to 
agree in general, but are not in good agreement. This is because the periodic 

terms (annual and Chandler) in pole coordinates ore not averaged out in the 

comparison duration (ten months). 

Finally, the transformation parameters between various reference frames 
fnherent in different techniques are calculated in terms of pole coordinates, 
using the following models: 

-Ay = -(yl - y^^) = -fix + «i cos • + a, sin • 

-Ax = -(x^ - = -fi 2 - «i sin • + aj cos • 

where fit, fit are the rotation angles about the first and second axes of the 
terrestrial coordinate system II; a,, a, of the corresponding inertial coordinate 
system 11; • is the Greenwich sidereal time. The values of fix, fij, a,, a, angles 
and their standard deviations with post-fit standard deviations of 

"observations" are listed in Table 3. 

It can be seen from Table 3 that the rotation angles are of the order of 10 
milliarcseconds, much smaller than expected. The interesting thing is that 
post-fit standard deviations (and thus the standard deviations of 

transformation parameters) for VLBI CSR39 are the smallest among them, while 
those associated with the BIH system are much larger. This is possibly an 
indication that VLBI and CSR39 pole coordinates are more consistent and 
probably more precis* . On the other hand, the BIH coordinates are seemingly 
questionable. Doppler coordinates are in between. 
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Table 1 Coefficients of Annual and Chandler Motions with Their Standard 
Deviations (units in (71001) 
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44 
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40 
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37 
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32 
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Table 2 


(units in 01001) 
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Table 3 (units in 01001) 
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